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Abstract 

A mean field argument is used to derive a master equation for systems simultaneously interacting 
with external fields and coupled environmental degrees of freedom. We prove that this master 
equation preserves positivity of the reduced density matrix. Solutions of the master equation are 
compared with exact solutions for a system consisting of three spins which is manipulated with a 
sequence of laser pulses while interacting with a spin-bath. Exact solutions appear to converge to 
the master equation result as the number of bath spins increases. 
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I. INTRODUCTION 



Decoherence constitutes a potentially serious problem for a wide range of proposed tech- 
nologies which exploit quantum phase interference. Examples include molecular electron- 
ics, laser control of chemical reactions, quantum computing and molecular motors. Since 
condensed phase implementation of such schemes is probably necessary for nontrivial appli- 
cations, some understanding of the influence of bath dynamics on system-bath interaction 
could prove useful for avoiding decoherence. Representation of condensed phase environ- 
ments with the uncoupled oscillator baths commonly assumed as a starting point in older 
theories[l|] is not generally valid. Indeed, recent experimental [2j and numericalpl] evidence 
supports the notion that atoms and molecules in condensed phases exhibit chaotic dynam- 
ics and hence strong bath self-interaction. In addition, some numerical simulations suggest 
that baths which have strong self-interaction cause much less decoherence Q than would be 
predicted by uncoupled oscillator models. Hence a theory of subsystem dynamics which 
can account for the effects of self-interacting baths could have important application in this 
area. Since lasers are often employed to initialize and manipulate such technologies the 
theory should also allow for interactions of the subsystem with external fields. 



We recently used a non-perturbative mean field approximation to derive a non-Mar 



Kovian 



positivity-preserving master equation for systems interacting with coupled bathsp, |6|, 
Tests of this equation against exact results for a spin interacting with a coupled spin-bath 
showed good agreement j^. Here we extend the theory by allowing the subsystem to interact 
with external time-varying electromagnetic fields. 

A natural starting point for any theory of subsystem dynamics is the exact projection 
operator approach introduced by Nakajima and Zwanzigj^. In the case of time-independent 
subsystem Hamiltonians this leads to a unique (up to the definition of the projection op- 
erator) and exact integro-differential master equation. Unfortunately, the memory kernel 
which weights the contributions of the previous states of the subsystem to its future cannot 



or a particular choice of projection 
8|. The extension to time-dependent 



be calculated in practice. Approximation of this kerne 
operator formed the basis of our previous workP, 0. 
subsystem Hamiltonians requires a similar derivation from first principles. This is outlined 
in section II. The essential step, as in our previous work,5^], is an appropriate mean field 
approximation for the interactions of the subsystem and environment. 
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The reduced density matrix of the subsystem should be positive semi-definite but violation 
of this property is a common failing of many master equations^]. We therefore prove in 
section III that the master equation obtained in section II preserves positivity. 

In section IV we introduce a model system which represents a set of three qubits of a 
quantum computer which are manipulated with a sequence of laser pulses. This three spin 
subsystem is allowed to interact with a bath of strongly coupled spins designed to model 
solid state environmental modes. Solutions are obtained for baths with varying numbers of 
spins. In section V we discuss the method used to obtain solutions of the master equation. 
The exact and master equation solutions are compared in section VI. We show that the 
exact results approach the mean field results as the number of bath modes increases. 



II. MEAN FIELD MASTER EQUATION 

Define a Nakajima-Zwanzigj^ projection operator P on the total (system plus bath) 
density xif) such that 

Px{t) = p{t)B, (1) 

where p{t) is the system density and B is the canonical bath density. Similarly, define 
Q = l-P. 

Consider a Hamiltonian 

H{t) =Hs + S{t) + H, + J2 S>.R>. (2) 

where S{t) involves only system operators and represents the effects of external electromag- 
netic fields. The operators Hg and represent the system while Hf, and represent the 
bath. Define a time-independent operator 

L = {l/h)[Hs + H, + Y^ S^R^, ■ ] (3) 

which is the Liouville operator in the absence of external fields, and a time- dependent 
operator 

m = m)mr] (4) 

for the lasers. Then using the facts that P and J^{t) commute and that PQ = it can be 
shown that 

dPx{t)/dt = -{i/n)[PLP + J^{t)] Px{t) - {i/h)PLQ Qx{t) (5) 
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dQx{t)ldt = -{i/h)[QLQ + Ht)] Qx{t) - {i/h)QLP Px{t). (6) 

Equation © should now be solved for Qx{t) and substituted into Eq. (jSJ. 

We have shown elsewhere js], 0] that P is non-Hermitian and that therefore QLQ is non- 
Hermitian with a complex spectrum. Let ujj and 7^ denote the real and imaginary parts of 



an eigenvalue of QLQ and let \(f)j) and be the associated right and left eigenvectors [11 1. 
We may expand Qx{t) = J2j Cj(t)\(j)j) in the complete eigenbasis and using orthonormality 
(i.e = Sj,k) it follows that ® can be rewritten in the form 

dC,it)/dt = - ^,)C,{t) - it/h)J2i^,\m\<l>k)Ckit) - i^/h)i^,\QLP Px{t). (7) 

k 

Now, the matrix elements (<l'j|jF(t)|0fc) should be dominated by the overlap of the bath 
part of the generalized eigenstates. For a sufficiently large bath the overlap represents an 
integration over a product of two essentially random functions. On this basis we should 
expect ($j|jF(t)|0jt) to vanish for j 7^ k. Secondly, note that since £{t) is Hermitian, 
has a real spectrum which is symmetric about zero, and hence the diagonal elements 
{(^j\T{t)\(j)j) are zero on average. The first approximation in our derivation is thus to neglect 
the matrix elements of J^{t) in Eq. (j?!). Assuming an initial state of the form x(0) = p(0)S 
we then obtain 

C^{t) = -{i/h) f dt'e^-'^^-^^^^'-''\^j\QLP Pxit') (8) 
which in turn then gives the desired solution 

Qx{t) = -{i/h) fdt'Y.e^-^-^-^^^^'-''%)i<l>,\QLP Px{t') (9) 
Jo ■ 

which can be substituted into equation 

Before making this substitution we require further approximation to eliminate the ex- 
plicit dependence on the unknown generalized eigenvalues and eigenvectors. Our previous 
approximation exploited the large number of bath degrees of freedom and the consequent 
randomness of the generalized eigenvectors. The large number of bath modes also implies a 
large spectral density of states for ujj and 7-,-. A large number of terms will thus contribute 
to the sum in (jH)) suggesting that perhaps the sum can be replaced by its average. That is, 
suppose that 



j j 



where W{t) = (cos cute"'''*) and the angle brackets denote an average over the generahzed 
spectral density. Here we have used the fact^, ^ that the spectral density is symmetric 
under u —u and the closure relation for the generalized eigenvectors. This mean field 
type approximation, which should be accurate for sufficiently large baths, implies that 

Qx{t) = -{i/h) f dt' Wit - t')QLP Pxit'). (11) 

Substituting (fTTj) into Eq. (0) and tracing over bath degrees of freedom it can then be 
shown i^i] that 

dp{t)/dt = -{t/h)[H, + S{t) + Y,R^S„p{t)] 

- ii/h')J2c,,, fdt' w{t-t'){[p{t')s,,s,] + [5., (12) 

where = TibiR^B} and C^j, = TTb{{R^ — R^){R^ — R^)B} denote canonical (i.e. B = 
^-Hh/kT ^r^^^^^-Ht/kT^^ averages and variances of bath operators. This is essentially the 
same master equation derived in Ref. 0] except the subsystem Hamiltonian is now time 
dependent. Note that W{t) plays the role of a memory function: it weights the contributions 
of previous states of the system. 

A careful treatment of the spectral properties of ^ = QLQ allows one to calculate the 
mean spectral density which in turn can be used to calculate W(t). Using this approach we 
have shown 0] that 

W{t) = [1 - ^{ptY + ^{ptr - ^{ptr + ^(pt)^]e-(''*)V« (13) 

where 



p = [{AA^) - {AA)]/J{AA^) (14) 



q = [{AA^) + {AA)]/^{AM). (15) 

The angle brackets here denote an average over the Liouville-Hilbert space i.e., for any 
operator F, 

m n 

= Ji^jy^^)T.T.(hj\F\hj) (16) 

i=l j=l 

where states denote a complete set. Simplified formulas for the real parameters (AA^) 
and {AA) are provided in Appendix A of Ref. [^. The memory function (fT!?|l is always 
positive and usually has a shape which is nearly gaussian. 

Note that the assumptions we have made in the above derivation mean that the master 
equation (|T^ is valid in the limit of large bath. 



III. PROOF OF POSITIVITY 



Define operators L(t) = {l/h)[Hs + £{t) + J^fiRf^S^,- ] and Ld such that L^p = 
C^,u{[pSu, S^] + [5^, 5^p]} and a function M{t) = r6{t) - W{t) where r = dt W{t). 
We may then write (|12|) in the form 

dp{t)/dt = -{iL{t) + TLd}p{t) + f'dt'M{t - t')Ldp{t'). (17) 

Jo 

Rather than attempt to prove positivity for p7|l we consider instead the related equation 
dpit, s)/dt = -{i{L{s) - id/ds) + rL4p(t, s) + f dt'M{t - t')Ldp{t', s) (18) 







in which we have introduced a new variable s which eliminates the explicit time dependence 
of L. [Note the similarity of this transformation to that employed in the (t,t') method jl^.] 
We will now show that (fTH|l preserves positivity of p(t,s), and since p(t) = p{t,s)\s=t, that 
(fTTjl preserves positivity of p{t). 

Consider that both L{s) and id/ds are Hermitian operators and that is of completely- 
positive-dynamical-semigroup form. It follows that the operator T> = —i{L{s) —id/ds) — 



rLd is the generator of a completely-positive-dynamical-semigroup jl3|]. Laplace transform- 
ing (fTHjl shows that p(t, s) = T(t, s)p(O) where the propagator T{t, s) is obtained by inverting 
the equation 

R{z,V + M{z)Ld) = {z-V-M{z)Ld)-' (19) 

= / dte-''T{t,s) (20) 
Jo 

for the resolvent operator R{z,V + M{z)Ld). Here M{z) is the Laplace transform of M{t). 

To invert Eq. ^lU^ it is convenient to first show that T{t, s) can be expressed in dynamical 
semigroup form. This means that there exists an operator A such that T(t, s) = e"^*. 
Equivalently, this means that the propagator can be written in terms of the resolvent via 



T(t, s) = lim 



Tl 

-R{n/t,A) 



(21) 



which implies that T{t,s) will preserve positivity if R{z,A) = {z — A)^^ is positive for 
large real z. Since the resolvent R{z, A) must be equivalent to the resolvent R{z, T> + 
M{z)Ld) (i.e. the solution p(t, s) is unique) it follows that T(t, s) will preserve positivity if 
R{z,V + M{z)Ld) is positive for large real z. This can be readily proved as we show below. 



Thus, the essential step is showing that the solutions of (fTHj) can be written in dynamical 
semigroup form. Substituting A = V + dt'M{t')LdS-t' into dp(t,s)/dt = Ap(t,s) and 
using the boundary condition p(t) = for t < one readily recovers Eq. (fTS|) . [Here 
6-t'f{t) = f{t — t') is the delay operator. A detailed derivation of A is presented in Ref. 
^.] Hence, the solutions of Eq. (|18|) can indeed be written in dynamical semigroup form 
and we need only show that our original resolvent is a positive operator. 
Clearly we may write 

R{z,V + M{z)Ld) = R{z,V) (l - M{z)LdR{z,V)y^ (22) 
= R{z,V)J2[m{z)L,R{z,V)] , (23) 

k=0 

and since both and R{z,V) = {z — V)^^ preserve positivity it follows that T[t,s) will 
preserve positivity if M{z) is positive for large real z. Now, since 

oo 

zt ; 



M{z) = / dte-'^M{t) (24) 

JO 

roo 

= T- dt e-''W{t) (25) 
Jo 

/■oo /"OO 

= / dt W{t) - / dt e-''W{t) (26) 
Jo Jo 

' dt {l-e-'')W{t) (27) 



it follows that M{z) is positive if W{t) is positive. In fact this is true for the memory 
function of equation (fT!?|) . Hence, Eq. (fTH|l preserves positivity of p(t,s). 
Finally, define p{t) = p{t, s)\s=t from which it follows that 

dp{t)/dt = dp{t, s)/dt\s=t + dpit, s)/ds\s=t (28) 

and inserting (fTHj) then gives our original equation (fTTj). Thus, positivity of p(t) is preserved 
by the master equation (fT^ . 



IV. SPIN-SPIN-BATH MODEL 



Our model system represents three qubits (i.e. two-level systems) which are manipulated 
with a sequence of laser pulses while the whole system interacts with an environment. Two 
of the qubits represent pairs of electronic states of impurities in a crystalline solid at low 
temperature. The third qubit represents two vibrational levels of an optical phonon mode 
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which is used as a means of transferring information from one impurity to the other via 
vibronic couphng induced by a laser. Other vibrational modes of the crystal (i.e. the 
environment) are represented by a number Ug of coupled spin-1/2 modes. Obviously the 

representation of phonon modes by spin-1/2 modes is valid only at low temperature. Here 
we set kT = .0067 eV which corresponds to liquid nitrogen temperature. 

FIG. 1: Without dissipation 
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We thus employ a time-dependent Hamiltonian of the form 
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ris t , ris-l ris 

j=l i=l j=i+l 

fc=0 

where we chose the gap between the electronic states of the impurities (labeled and 2) to 
be hujeg = 3 eV. The excitation energy of the optical phonon (labeled 1) was set to hujp = .2 
eV. A large system-bath coupling of hXo = .0075 eV was chosen so that decoherence could 
be observed during the roughly 160 fs of time evolution. An intra-bath coupling of hX = .03 
eV was chosen which is roughly representative of diamond. A small anharmonic term with 
hf3 = .0001 eV was included. Bath frequencies were sampled from a Debye distribution with 
a cutoff at hujc = .05 eV. 

The laser interactions move an initial excitation from the first impurity to the optical 
phonon and then from the optical phonon to the second impurity. The process is then 
reversed. Overall we repeat this sequence three times. The parameters of the lasers are 
^^laser = ^eg — (Igss than the diamond band gap of 5.4 eV), ha = .325 eV, b = .325a^. 
Finally, the first pulse sequence times are ti = 10 h/eV, t2 = 30 h/eV, = 50 h/eV and 
^4 = 70 ^/eV with multiples delayed by r = 80 h/eV. 

We calculated the reduced density of the three qubit system via the formula 

pit) = Y.Pm Trb{|^„(t))(^„(t)|} (30) 

m=l 

where 

'^eig 

Pm = exp{-em/kT}/ ^ exp{-ei/kT}, (31) 

1=1 

em and |m) are bath energies and eigenvectors (i.e. of terms 5 and 6 of Eq. ()29|)). and 
kT = .0067 eV is the (liquid nitrogen) temperature in units of energy. The notation 
T^T^b{\i'm(t)) {ipm(t)\} iudicatcs a trace of the full density \ipm(t)) {4'rn(t)\ over the environ- 
mental degrees of freedom. The states \ipm{t)) are evolved via the Schrodinger equation 
from initial states 

IV^^(O)) = |100)® |m) (32) 
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under Hamiltonian (j^ . The basis of eigenstates of the Oz operators was used to represent 
all states. The state |100) in (jH^ refers to the system and means that the 0-spin was initially 
excited while the 1-spin and 2-spin were in their ground states. 

The ARPACK linear algebra software [l^ was used to calculate the lowest ngjg = 10 ener- 
gies and eigenvectors of the isolated environment. The temperature was chosen such that no 
states with quantum number m higher than n^ig are populated at equilibrium. The numer- 
ical solutions of the Schrodinger ordinary differential equations for were calculated 
using an eighth order Runge-Kutta routine^^. Operations of the Hamiltonian ()29|) on the 
statevector were calculated via repeated application of Pauli matrix multiplication routines. 
For example 

Oo, Jl, • • • , Ji, • • • , Jn.+2ki'^|V^) = Oo, jl, • • • , Ji, • • • ,in,+2|'0) (33) 

for all sets of j/ = 0, 1, / = 0, 1, . . . , + 2 and where jj = 1 if jj = and jj = if jj = 1. 
Thus, an operation of a^^ simply rearranges the components of States of the basis 
can be represented by integers j = jo + Ji2 + + . . . + jj2* + . . . + and since 

integers are represented in binary form on a computer, the mapping j j' = jo + ji2 + 
j22^ -|- . . . + jj2* + . . . + jn^+22"''"'"^ under aj^^ can be calculated very simply using Fortran 
binary-operation system functions. Operations for a^^^ and o"^*-* are also straightforward. 

Finally, from the reduced density (jHn|) of the system we calculated the reduced densities 
of the two qubits and of the optical phonon mode by tracing out the remaining unwanted 
degrees of freedom. For example, for qubit we calculated 

pW(t) = TriTr2{p(t)} (34) 

while for qubit 2 and for phonon mode 1 

p(2)(t) = TroTri{p(t)} (35) 
p(i)(t) = TroTr2{p(t)}. (36) 

Thus, each component of the system is represented by a 2x2 matrix which makes it easier 
to display the solutions and compare them with solutions of the master equation. 

In order to show convergence to the master equation results we considered a range of 
values of rig. Specifically, we report results for Ug = 10, 12, 14 and 16. 
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For reference we plot the solutions for the subsystem in the absence of dissipation in 
Fig. 1. Figure 1(a) shows the occupation probabilities for the ground state (solid curve) 
and excited state (dashed) of the first impurity plotted against time in units of h/eV= .66 
fs. The real and imaginary parts of p^^i are identically zero. Similar quantities for the 
optical phonon and the second impurity are plotted in 1(b) and 1(c), respectively. Again 
the off-diagonal elements are zero as a consequence of our choice of initial state. 



V. NUMERICAL SOLUTION OF MASTER EQUATION 

We recently developed a numerical technique for solving integro-differential equations 
like fjl2j) . The accuracy of the method has been established for both generalized Langevin 
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equations and master equations by comparison with exact solutions[17|. Basically the 
method works by converting integro-differential equations to ordinary differential equations 
which are solved by standard methods. 

We implement the method as follows. Define a space-like time variable u and a smoothed 
density operator 

X(t, u) = fiu) f dt' Wit -t' + u)p{t'), (37) 
Jo 

where f{u) is a damping function such that /(O) = 1. Direct substitution shows that p(t) 
and x(^; u) satisfy ordinary differential equations 

dp{t)/dt = -{i/h)[H, + £{t) + Y.R,S^,p{t)] 

- {i/h^)J2c^A[x{t,o)s,,s,] + [S,,S,x{t,0)]}, (38) 

fJ,,U 

dxit, u)/dt = f{u)W{u)p{t) + dxit, u)/du - U\u)/f{u)) x{t, u) (39) 



(where f'{u) = df{u)/du) which are then solved by representing m on a grid. 
More specifically, the equations for Hamiltonian (|29p are 

dp{t)/dt = -%\ujj2af + uj^llaf^ + uj^gjlaf + Ao(a(°) + ai^^ + af^)Y., 

fc=0 

+ a(^Vfe-''(*-*-'^-)^+aWa«e-^(*-*-^^)^}cos^,.,,.t, p{t)\ 



- C{{af + + a(^))^x(t, 0) + Xit, 0)(a(°) + a^J^ + af^f 
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FIG. 2: Impurity 1 with dissipation 
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dx{t, u)/dt 



W{u)p{t) + dx{t, u)/du + 2gu x{t, u) 



(40) 
(41) 



where = Tn{J:,B}, = E"=i and C = A^Tr^lS^S}. 

Equations (jiUll and (PT|) were restricted to a grid of points Uj = (n + / — j)At with 
, = i = »*M38«) whe.e A* = .1 R/eV . the time-step employed m the 

dynamics. Following Ref. |17[ a damping function f{u) = e ^" with (7 = ll/[{n — l)At]'^ was 
used. Converged results were obtained for n = 100 grid points. We chose W{u) = W{\u\) for 
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negative values of u. A discrete- variable jiol matrix representation was employed to calculate 
the partial derivative with respect to u in Eq. pTj) . 

FIG. 3: Optical phonon with dissipation 
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Fmally, the ordmary d.ffe^nt.al equations m and m were mtegrated usmg an e.ghth 
order Runge-Kutta routine jl5|. 

The parameters of the master equation and C, and the parameters of the memory 
function (jl3|) . were calculated using the exact energies and eigenvectors of the bath Hamil- 
tonian computed in Section III and formulas reported in Ref. j^. The parameters of the 
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master equation converge rapidly with the number of bath spins since they are all aver- 
age quantities. Thus, master equation solutions for Us — 12 cannot be distinguished from 
solutions with = 14 or Ug — 16. 

VI. RESULTS 

Figures 2-4 show the occupation probabilities and real and imaginary parts of off-diagonal 
density elements for each of the three qubits, for various values of the number of bath modes 
function of time. Curves for Ug — 10 (dashed), Ug — 12 (short-dashed), Ug — 14 
(dotted), and Ug — 16 (dot-dashed) are shown in each figure. For comparison we also show 
the solutions of the master equation (solid curve). [Note that the occupation probabilities 
for the master equation are positive, in agreement with our proof in section III.] Time t is 
in units of /i/eV~ .66 fs. 

Each sequence of four laser pulses can be viewed as moving an excitation from the first 
impurity to the optical phonon and from the phonon to the second impurity, then back to the 
phonon and finally back to the first impurity. The sequence is repeated three times for a total 
of twelve laser pulses. The general idea is to roughly simulate the sort of manipulations that 
would be employed in a quantum computer. Because calculation of the system dynamics 
in the presence of the bath spins is very expensive we have chosen a strong system-bath 
coupling and short pulse width so that decoherence is manifested over the relatively short 
time span of 160 fs. 

In accord with the initial conditions and pulse sequence Fig. 2 shows an excitation on 
impurity 1 {pfi = 1, Pqq^ = 0), which relaxes to its ground state {pfi = 0, Pqq^ = 1) and 
then is rc-excitcd. This is repeated three times. Figure 3 shows the phonon mode initially 
in its ground state {p^ii = 0, Pqq^ = 1). The phonon is then excited (p^^"* = 1, p^^Q = 0) 
and then returned to its ground state. Again, this is repeated three times. The second 
impurity, shown in Fig. 4, is initially in its ground state (p^i — 0, poo^ = 1)- It is then 
excited (p^i — 1, p^Q — 0) and then returned to its ground state. The real and imaginary 
parts of the off-diagonal elements of the three qubits show small oscillations throughout the 
manipulations as a result of decoherence. These oscillations are much faster for the two 
impurities and so we show only the last forty time units. Compare Fig. 1 with Figs. 2-4 
and note the obvious effects of decoherence and dissipation present in the dynamics of all 
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density components. Return to the the initial state is imperfect because of decoherence. 
Note also that while the off-diagonal elements are zero in the absence of dissipation, here 
they show small oscillations. 



FIG. 4: Impurity 2 with dissipation 
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It is clear that the agreement between the exact calculations and the master equation 
improves dramatically with increasing numbers of bath spins. Decoherence effects are gen- 
erally much stronger for small baths but these decrease as the bath gets larger. For small 
numbers of bath spins, the return to the initial state after a pulse sequence is less perfect 



15 



than that predicted by the master equation. As the number of bath spins increases this 
discrepancy is incrementally reduced. Deviations of the real and imaginary parts of the 
off-diagonal elements from master equation predictions also decline as the number of bath 
spins increases. 

The results for 16 bath spins show close but still imperfect agreement with the predictions 
of the master equation. It is obviously of interest whether the exact numerical solutions 
converge precisely to the master equation results in the limit that the bath is very large. 
Unfortunately, we were unable to perform calculations for larger baths and so this will have 
to remain an open question. 

VII. SUMMARY 

We have used a mean field approximation to derive a master equation suitable for time- 
dependent subsystem Hamiltonians and self-interacting baths. After proving that the master 
equation preserves positivity of the reduced density matrix we compared its solutions to those 
of a model system. We found that exact numerical solutions for the model converged toward 
those of the master equation as the number of bath modes was increased. This supports our 
expectation that the approximate master equation will become increasing accurate as the 
bath size approaches the thermodynamic limit. 

We have recently developed an exact method for decomposing the quantum N-body 
vibronic dynamics problem (for pairwise interaction) into N stochastic 1-body problems jl8|. 
That is, we can now exactly solve the dynamics of pairwise interacting distinguishable spins 
and vibrations. This should allow us to obtain exact solutions for more realistic models and 
for larger baths. We hope to soon test the predictions of the mean field master equation 
against exact solutions for these more realistic models. 

The authors gratefully acknowledge the support of the Natural Sciences and Engineering 
Research Council of Canada. 
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